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ABSTRACT 

An explicit finite difference real time iteration scheme is developed to study harmonic sound propagation in 
aircraft engine nacelles. To reduce storage requirements for future large three-dimensional problems, the time 
dependent potential form of the acoustic wave equation is used. To insure that the finite difference scheme is both 
explicit and stable for a harmonic monochromatic sound field, a parabolic (in time) approximation is introduced to 
reduce the order of the governing equation. The analysis begins with a harmonic sound source radiating into a quies- 
cent duct. This fully explicit iteration method then calculates stepwise in time to obtain the “steady state” harmonic 
solutions of the acoustic field. For stability, applications of conventional impedance boundary conditions requires 
coupling to explicit hyperbolic difference equations at the boundary. 

The introduction of the time parameter eliminates the large matrix storage requirements normally associated 
with frequency domain solutions, and time marching attains the steady-state quickly enough to make the method 
favorable when compared to frequency domain methods. For validation, this transient-frequency domain method is 
applied to sound propagation in a two-dimensional hard wall duct with plug flow. 


INTRODUCTION 

Both steady-state (frequency domain) and transient (time domain) finite difference and finite element tech- 
niques have been developed to study sound propagation in aircraft nacelles. To date, the numerical solutions have 
generally been limited to moderate frequency sound and mean flow Mach numbers in two-dimensional axisym- 
metric nacelles. Wavelength resolution problems have prevented a broader range of applications of the numerical 
methods. A fine grid is required to resolve the short wavelengths associated with high frequency sound propagation 
with high inlet Mach numbers. Thus, application of numerical techniques to high frequency sound propagation in 
three-dimensional engine nacelles has yet to be attempted. 

To extend numerical analysis to higher frequencies and inlet flow Mach numbers, as well as three-dimensional 
geometries Baumeister and Kreider (ref. 10), suggest using a pseudo-time-frequency transformation to the acoustic 
potential equations. Their method eliminates the large matrix storage requirements of the steady state (Fourier trans- 
forms) techniques in the frequency domain but still allows the use of conventional impedance conditions. Most 
importantly, their two step time marching formulation is fully explicit under flow conditions. They also suggested 
using an almost highly structured grid to reducing computer storage and run times. Using the same grid system and 
governing acoustic potential equations, the goal of the present paper is to develop a stable real time iteration scheme. 
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The paper begins with a brief description of the governing equations. A detailed development of the governing 
equations has been given earlier (ref. 10); it is therefore described only briefly here. Next, assuming a harmonic 
monochromatic sound field, a parabolic (in time) approximation is introduced to reduce the order of the governing 
equation. The bulk of the paper describes the development of a stable, explicit finite difference scheme. The scheme 
is iterated in time to converge to the steady-state solution associated with a Fourier transform solution. 

NOMENCLATURE 

c steady speed of sound, C # /C^, equation (2) 

D # dimensional duct height or diameter 
D duct height, D = 1 
d parameter, equation (25) 

F source amplitude at duct entrance, equation (17) 
f* dimensional frequency 
f dimensionless frequency, f^D^C* 
g parameter, equation (27) 
h parameter, equation (26) 

i -\/— T 

L length of duct, L # /D # , figure 1 
M f Mach number at duct entrance 
n unit outward normal 
P^ dimensional pressure 
P dimensionless fluid pressure, PfypJC^f 2 
P' acoustic pressure fluctuation, equation (6) 
t dimensionless time, f% # 

At time step 

x dimensionless axial coordinate, x # /D # 

Ax axial grid spacing 

y dimensionless transverse coordinate, y # /D # 

Ay transverse grid spacing 
y ratio of specific heats 
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p dimensionless fluid density, p # /p J 
p steady fluid density 

0 dimensionless time dependent flow potential, (J) # /C^D # 
(j> steady mean flow potential 

(f)' transient acoustic potential 
\\f spatial potential, equation (10) 
co dimensionless frequency, 2nf 
V D # V # 

Subscripts 

1 axial index, see figure 1 

j transverse index, see figure 1 
o ambient of reference condition 
Superscripts 

# dimensional quantity 
k time step 


PROBLEM STATEMENT 

The problem under consideration here is the steady-state propagation of sound, represented by the perturbation 
acoustic potential, through a two-dimensional rectangular duct. The source, noise emanating from fan blades in a jet 
engine inlet nozzle, is represented by specifying the pressure distribution at the fan face. The goal of the paper is to 
develop a stable, explicit finite difference scheme that incorporates the far field impedance condition applied at the 
duct exit and rigid body boundary conditions on the duct walls. The method is designed with the intention of extend- 
ing the current two-dimensional duct formulation to general three-dimensional nacelle design problems with a vari- 
ety of possible boundary conditions in the near and far fields. 


GOVERNING EQUATIONS 

Acoustic propagation in inlet nacelles can be reasonably modeled by an inviscid approximation. For single 
mode JT15D engine data, a previous finite element study (ref. 8) employing the potential formulation in the fre- 
quency domain showed good agreement with experimental data — in the far field radiation pattern as well as sup- 
pressor attenuation. Due to this success, the problem under consideration here is formulated in terms of an acoustic 
potential. 

For inviscid, nonheat conducting and irrotational flow, the linearized acoustic equation for two-dimensional 
potential flow can be written in dimensionless form as (ref. 10, eq. (1)) 
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0 = fVtt - (c 2 - ®x )<t> X x - if - fyjVyy + 2® x ® y <|)' Iy + 2fO x <(> xt + 2f4> y <|> yt + 2(® x ® xx + <t> y O xy )«t>' x 

4 2|O x O xy + ®y®yy j^y — (y ~ 0(^t ^x^x 4 ^y^y ji^xx 4 ®yy) (1) 

where 

5 2 = l-|(Y-l)(®x+®y) (2) 

The symbol <j> represents steady mean flow potential while cj)’ represent the time dependent acoustic potential. The 
speed of propagation of a disturbance is represented by c and the frequency of an acoustic source by f. The sub- 
scripts indicate partial differentiation with respect to the subscripted variables. 

The conventional normalization factors used to develop these nondimensional equations are given in the 
NOMENCLATURE. However, the normalization of time deserves some special comment. A common choice for 
normalizing time is t = C^t # /D # . The superscript # designates a dimensional quantity while the subscript o indicates 
an arbitrary reference value. With this choice, the dimensionless frequency f would not appear in equations (1) or 
(2). However, in this paper, the dimensional frequency f* of the forcing acoustic signal was chosen to normalize 
time, so that t = f^t # . As a result, the time t indicates the number of complete acoustic cycles that have occurred since 
the start of the solution process. This is advantageous because the total time of the numerical calculation can gener- 
ally be set independently of the frequency of the acoustic signal. 

For simplicity, the current paper will focus only on plug flow. For this case, the mean flow terms in equation (1) 
become 


O y = <D xy = O xx = 0 O x =M f (3) 

Substituting equation (3) into equation (1) yields 

0 - f : V tt - (c 2 - Mf )<t>' xx - c 2 <t> yy + 2fM f <j>' xt (4) 

and 

c 2 =1_I( y -1)m? (5) 

The relationship between the acoustic pressure and potential can also be expressed as (ref. 10, eq. (7)) 

ip'(x,y,t) = -f<|)' t -M t (|> x (x = 0, ® y =o) (6) 

Equation (4) is the basic equation used to establish a finite difference formulation for sound propagation in the rect- 
angular duct shown in figure 1. However, care must be taken when discretizing derivative terms to insure that the 
resulting scheme is both stable and explicit. Note that it is easy to develop a stable implicit method, but this would 
yield a matrix formulation that is no better than a frequency domain approach. It turns out that the proper treatment 
of the mixed time and space derivative terms (which appear on the second line of eq. (1)) is critical for maintaining 
stability. 

The implicit formulation is obtained by approximating the mixed partials by (ref. 6, eq. (16), and ref. 1, p. 884, 
eq. (25.3.27)) 


<l>xt = 


x k+1 . i' k i /k' k i a k — 1 ni k a' k+1 ^ k — 1 

foi.j + Qi-lJ + foi+l.j + QiJ ~ ^foi,j ~ foi— l,j ~foi+l,j 

2AxAt 


( 7 ) 
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( 8 ) 


^yt ~ 


j.' k+1 , j.' k- 

_ + 4>i,j 


4 + + <i>'i,j-l - 24>if - 4>'i j+i 1 - Vi,] 1 -! 1 

2AtAy 


where i and j denote the space indices for the nodal system shown in figure 1, k is the time index defined by 

t k+1 = t k + At (9) 

and Ax, Ay, and At are the space and time mesh spacings, respectively. 

Reference 6 shows that equation (7) can be used in an explicit fashion for one-dimensional plug flow problems 
in a duct. However, if the flow is not one-dimensional or if the region exterior to the duct is included, the scheme 
must be implicit. This approach is inappropriate for general three-dimensional problems. This problem can be cir- 
cumvented by modifying the governing equation. The details follow in the next section. 


PARABOLIC APPROXIMATION 

There are several ways to develop a frequency domain formulation for the general two-dimensional acoustic 
wave equation (1) or the plug flow simplification equation (4). The Fourier Transform can be applied if the potential 
has a multifrequency content. In the monochromatic case, this is equivalent to assuming that 

(}) , (x,y,t) = \|/(x,y)e _ia) 1 = \|/(x, y)e -l27Ct (10) 


which, in the case of plug flow (from eq. (4)), yields 

0 = (c 2 -Mf )\|/ xx + c 2 \|/ yy + co 2 \|/ + i2coM f \|/ x (11) 

This equation would be solved numerically using a linear system of equations. However, the associated matrix 
is not positive definite, which can lead to numerical difficulties, and which preclude the use of iterative techniques. 
Therefore, it is desirable to develop an explicit finite difference scheme to avoid the use of matrices. In time- 
dependent form, equations (1) or (4) cannot easily be discretized in such a way that the resulting finite difference 
scheme is both stable and explicit in the presence off low (it is possible to obtain reasonable results in the no-flow 
case). 

Reference 10 resolves these difficulties by a transformation to a transient-frequency domain. Herein, however, 
the resolution of these difficulties is achieved by applying the transformation equation (10) only once to equa- 
tion (4). Employing equation (10), the time derivatives in equation (4) can replaced by the following relationships: 

4>t = -i27E\|/e" i27Ct = -i2my (12) 

<t>xt = ^-(-i27C\|/e _l25tt ) = -i2w x e“ l23tt = -i27C<j> x (13) 

Under this transformation, the plug flow equation (4) becomes, 

— ifco(j) t = |c 2 — M 2 j(|) xx + c 2 (|) yy + i2coMf <|) x (14) 

The mixed derivative problem term which prevented an explicit finite difference representation of equation (4) 
has been eliminated from equation (4). An explicit finite difference solution can now be formulated to solve 
equation (14). 
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The interest in this paper is in steady state harmonic solutions to the wave equation as represented by the vari- 
able V|/. Therefore, after sufficient time during the transient solution of equation (14), the ((>' variable can be related to 
the steady state Fourier transformed solution, using equation (10). 


Y(x, y) = 


4>'(x,y,t) 

e - i27It 


(15) 


INITIAL AND BOUNDARY CONDITIONS 

The duct is assumed to be quiescent at time 0, so that the initial condition is 

<j>' (x, y, 0) = 0 (16) 

As the equation is iterated in time, the solution builds up to the steady state harmonic solution. 

At the duct entrance, (x = 0), the potential is given by 

(0, y, t) = F(y)e~' 2m (17) 


If the pressure at x = 0 is specified as the boundary condition, then the potential is related to the pressure directly 
through equation (6). 

The hard wall condition is 


V<t>* • n = 0 


(18) 


where n is the unit outward normal. 

To simulate a nonreflective boundary at the duct exit in figure 1, the difference equation at the end of the com- 
putational domain can be expressed in terms of an exit impedance. In tern, this can be used to develop an exit gradi- 
ent condition to simulate a nonreflecting grid termination. For plane wave propagation in the examples to follow, 
Baumeister and Kreider (ref. 10) have shown this gradient to be of the form: 


<f>x = 


1(0 


c + lVL 


-v 


(19) 


FINITE DIFFERENCE EQUATIONS 

The finite difference approximations determine the potential at the spatial grid points at discrete time steps 
t k = kAt. Starting from the known initial conditions at t = 0 and the boundary conditions, the algorithm marches the 
solution out to later times. 

Away from the duct boundaries, as shown by the cell in figure 1, each partial derivative in equation (14) can be 
expressed as follows: 


i ' k+1 i ' k-1 

1 2At 


<i> xx = 


Ax 2 


( 20 ) 


( 21 ) 
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( 22 ) 


“ Ay 2 


<l>x = 


2Ax 


(23) 


Substituting these expressions into equation (14) yields 




2d 


2c 2 ^ 


Ax 2 Ay 2 


+ < t’i+ k l,j 


' d " + - g 


Ax 2 2 Ax 


+ ^i-lj 


"d 2 


g 


Ax 2 Ax 


+ ^i.j+1 




Ay 


, a k 

+ 4>i, j— 1 


J 


( t\ 

vAy 2 , 


H-V.f-'gl (24) 


where 


d 


2 



(25) 


and 


h = icof 


(26) 


and 


g = 2icoM f (27) 

Equation (24) is an explicit two step scheme. At t = 0, field values at t k_1 are assumed zero because the initial field is 
quiescent. 

The expressions for the difference equations at the boundaries are complicated somewhat by the impedance 
conditions. However, a simple integration procedure resolves this problem. Baumeister (refs. 5 and 6) gives precise 
details for generating the time difference equations at the boundaries. In fact, because of problems with applying the 
exit termination condition, the actual hyperbolic exit difference equations used in these references will be applied. 
This will be fully discussed later. 


STABILITY 

A von Neumann stability analysis (ref. 14) indicates that the method is conditionally stable, subject to the con- 
servative condition 


4 


fc] 

2“ 

, 2|M f | 

cof 

VAxJ i 

UyJ 


fAx 


( 28 ) 
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In a typical application, go, f, and M f are set by the engine operating conditions. Next, the grid spacing parameters 
Ax and Ay are set to accurately resolve the estimated spatial harmonic variation of the acoustic field. Finally, At is 
chosen to satisfy equation (28). 

In the von Neumann analysis, conditional stability means that the amplification factor, which describes how 
errors propagate from one time step to the next, has magnitude one. Thus, when inequality (28) is satisfied, errors 
are not magnified or diminished in magnitude. This is a desirable property, since the numerical formulation can not 
distinguish between an error and a small acoustic mode. 

The von Neumann stability analysis does not take into account boundary conditions. For stability, gradient 
boundary conditions require the use of smaller At than predicted by equation (28). 


NUMERICAL EXAMPLES 

In the three examples that follow, the parabolic transient results are compared to the exact results of the steady 
Fourier transformed solutions. The exact Fourier transformed solution for the problem considered is given by 
(ref. 10, eq. (35)) 


ico 

\K(x) = e" +MfX 


(29) 


The following problem is considered: a plane wave propagates from the left into a quiescent duct of length one, 
and the acoustic potential field is to be computed in the duct. Note that, boundary conditions can introduce instabili- 
ties (refs. 7 and 11) into otherwise stable finite difference schemes. Therefore, it is important to test the proposed 
method for convergence in time to the steady state solution in the absence of the exit boundary condition (eq. (19)), 
and to test independently the effect of the exit boundary condition itself on the solution. 


Semi-Infinite Duct 

In this example, the computational boundary is set far enough away from the true boundary x = 1 that any arti- 
facts arising from imperfections in the exit boundary condition do not affect the solution in [0,1]. The numerical 
solution propagates one node per time step, so setting the boundary at x = 50 with step Ax = 0.05 provides a suffi- 
cient number of time steps to gauge the convergence of the method before any artifacts might reflect back from the 
computational boundary. 

The numerical and exact results are compared in figure 2, for no flow (fig. 2(a)) and for Mach number 
M f = -0.5 (fig. 2(b)). The frequency is normalized to f = 1. Both cases show excellent agreement. The total calcula- 
tion time was tp - 5.0. 


Finite Duct L = 1 

In this example, the computational boundary is moved up to the true boundary x = 1 to examine the effect of 
the exit boundary condition (eq. (19)). The frequency is normalized to 1. Three cases were considered — no flow 
(M f = 0) and Mach number M f = ±0.2. 

The formulation using the parabolic wave equation, equation (14), was found to be unstable. Apparently, the 
parabolic wave equation is not compatible with impedance boundary conditions represented by equation (19). Fortu- 
nately, the exit termination difference equation could be written in an explicit fashion using the hyperbolic equation, 
equation (4). The difference equation for this exit condition is given by Baumeister (ref. 5). 

The results are shown in figures 3(a) to (c), respectively using a parabolic difference formulation in the interior 
of the duct and a hyperbolic formulation at the exit. In figure 3, the direction of the arrow indicates the direction of 
propagation of the acoustic wave. The numerical results again match well with the exact results. Notice also that the 
time step has been decreased here, which tends to increase the execution time; however, the computational domain 
is smaller, which tends to decrease the execution time. The total calculation time in this example was tp = 4.0. 



It is clear that the exit boundary condition does have a slight degrading effect on the solution particularly with 
flow. In fact, for Mach number greater than 0.2 the errors are unacceptable. Therefore, the parabolic solution should 
only be employed in problems where the exit Mach number is zero or reasonably low. Consequently, in most prob- 
lems, the transient-frequency approach develop by Baumeister and Kreider (ref. 10) is preferred. 


CONVERGENCE RATE 

In this example, the convergence rate is studied for the region x = 0 to x = 1 using the semi-infinite duct. The 
results are shown in figure 4 for the magnitude of the potential. As seen in this figure, the numerical solution quickly 
and accurately converges to the exact steady state solution. As seen in figure 4, the parabolic solutions converges to 
the Fourier transformed results after a time of t = 1 has elapsed. This is about twice as fast as the transient-frequency 
approach presented by Baumeister and Kreider (ref. 10). 


SOLUTION METHODS 

With the parabolic transient approach developed in this paper, four different solution techniques are now avail- 
able to solve the hyperbolic wave equation that describes acoustic propagation in ducts and jet engine nacelles with 
a monochromatic source. This is illustrated in figure 5 for the zero mean flow case. 

The Fourier transform of the wave equation was the first numerical approach used to study sound propagation 
in jet engine ducts (refs. 2 and 4). This steady state approach is outlined on the centerline of figure 5 (third column 
from left). The governing hyperbolic wave equation is transformed to the elliptic Helmholtz “wave” equation. Finite 
difference (FD) and finite element (FE) numerical formulations have been employed to solve this equation. After 
applications of the boundary conditions (fig. 5; [BC]), the associated finite difference or finite element global matrix 
is solved for the velocity potential (or pressure). Because the matrix form of the Helmholtz partial differential equa- 
tion is not positive definite, matrix elimination solutions are generally employed. This require extensive storage. 
Conveniently, the steady state approach allows the direct calculation of the sound pressure levels. Currently, this is 
the most popular approach for solving acoustic propagation in ducts. 

In the inlet to a turbojet engine, the dimensionless frequencies f can be on the order of 30 to 50 for the higher 
harmonics of the blade passing frequency. The storage requirements and associated computer run times for these 
high frequencies makes computations expensive and even impossible. To make the numerical solutions more cost 
effective, grid saving approximations to the governing Helmholtz wave equation have been used (fig. 5; [Approx. 
Spatial]), Baumeister (ref. 3) employed the wave envelope theory while Hardin and Tappert (ref. 13) developed a 
similar approach for underwater sound propagation with the addition of a parabolic (space) approximation. Candel 
(ref. 12) presents an extensive discussion of the contemporary research in this area and a detailed development of the 
parabolic (spatial) equation method (PEM). 

The transient solution to the wave equation was the second numerical approach used to study propagation in jet 
engine ducts, which is shown in the second column from the left in figure 5. To eliminate the matrix storage require- 
ments, Baumeister developed time dependent finite difference numerical solutions for noise propagation in a two- 
dimensional duct (ref. 5). 

Sound is introduced as a boundary condition at the duct entrance. The initial conditions generally assume a qui- 
escent duct. Finite difference (FD) approximations to the hyperbolic wave equation are then solved by an iteration 
process. The calculations are run until the initial transient dies out and steady harmonic oscillations are established. 
Finally, the transient variable <t>’ is transformed into the steady state variable \j/ associated with the solution of the 
Helmholtz equation. As with the steady state approach, approximate spatial solutions reduce computer storage and 
run times (ref. 9). 

The third option, the transient-frequency technique (ref. 10), is illustrated in the fourth column from the left of 
figure 5. This fully explicit iteration method eliminates the large matrix storage requirements of steady state tech- 
niques and allows the use of conventional impedance conditions. As time increases, the iteration process directly 
computes the steady state variable \| /. 

The fourth option, as discussed in this paper, involves a parabolic (time) approximation to the wave equation. 
The transient finite difference solution of this equation can be conveniently expressed in explicit form with mean 
flow. 
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CONCLUDING REMARKS 


A parabolic transient numerical solution of the potential acoustic equations has been developed. The potential 
form of the governing equations has been employed to reduce the number of dependent variables and their associ- 
ated storage requirements. The method eliminates the large matrix storage requirements of steady state techniques in 
the frequency domain and the formulation is fully explicit under flow conditions. The field is iterated in time from 
an initial value of 0 to attain the steady state solution. In each example provided, the numerical solution quickly and 
accurately converges to the exact steady state solution. 

Application of impedance boundary conditions to this parabolic formulation causes an instability. To circum- 
vent this problem, the termination boundary conditions are developed in terms of the hyperbolic acoustic equations. 
However, with high Mach numbers this termination procedure leads to large errors. Thus, the method should only be 
applied to problems with a zero mean flow far field boundary condition. 
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Figure 3.™ Analytical and numerical potential profiles along wall for plane wave propagating in a hard 
wall duct of unit length and a non-reflecting exit condition (f - 1). (a) M f = 0 (Ax = 0.05, At = 0.001 t T = 4.0). 
(b) M f = +0.2 (Ax = 0.05, At = 0.001 t T = 4.0). (c) M f = -0.2 (Ax = 0.05, At = 0.001 t T = 4.0). 
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Figure 5. — Alternate finite difference/element methods in solving wave equation with monochromatic source. 
(1C - initial condition, BC - boundary condition). 
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